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I . INTRODUCTION 

The Penn State Finite Difference Time Domain Electromagnetic 
Scattering Code Versions TEA and TMA are two dimensional 
numerical electromagnetic scattering codes based upon the Finite 
Difference Time Domain Technique (FDTD) first proposed by Yee [1] 
in 1966. The supplied version of the codes are two versions of 
our current two dimensional FDTD code set. This manual provides 
a description of the codes and corresponding results for the 
default scattering problem. The manual is organized into eleven 
sections: introduction, Version TEA and TMA code capabilities, a 

brief description of the default scattering geometry, a brief 
description of each subroutine, a description of the include 
files ( TEACOM . FOR TMACOM. FOR) , a section briefly discussing 
scattering width computations, a section discussing the 
scattering results, a sample problem setup section, a new problem 
checklist, references and figure titles. 

II. VERSION TEA AND TMA CODE CAPABILITIES 

The Penn State University FDTD Electromagnetic Scattering 
Code Versions TEA and TMA have identical capabilities except the 
TEA code has electric field perpendicular to the z axis and the 
TMA code has electric field parallel to the z axis. Each code 
has the following capabilities: 

1) Ability to model lossy dielectric and perfectly conducting 
scatterers . 

2) First and second order outer radiation boundary condition 
(ORBC) operating on the electric fields for dielectric or 
perfectly conducting scatterers. 

3) Near to far zone transformation capability to obtain far zone 
scattered fields. 

4) Gaussian and smooth cosine incident waveforms with arbitrary 
incidence angles. 

5) Near zone field, current or power sampling capability. 

6) Companion codes for computing scattering width. 

III. DEFAULT SCATTERING GEOMETRY 

The codes as delivered are set up to calculate the far zone 
backscatter fields for an infinite, 0.25 m radius, perfectly 
conducting cylinder. The problem space size is 201 by 201 cells 
in the x and y directions, the cells are 1 cm squares, and the 
incident waveform is a Gaussian pulse with incidence angle of 
0=180 degrees. The output data files are included as a reference 
along with codes ( SWTEA . FOR , SWTMA. FOR) for computing the 
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frequency domain scattering width using these output data files. 
The ORBC is the second order absorbing boundary condition set 
forth by Mur [2]. 

IV. SUBROUTINE DESCRIPTION 

In the description for each subroutine, an asterisk (*) will 
be placed by the subroutine name if that particular subroutine is 
normally modified when defining a scattering problem. Also, each 
subroutine will be denoted if it is applicable to the TE code 
only, the TM code only, or to both codes. 

MAIN ROUTINE (TE, TM) 

The main routine in the program contains the calls for all 
necessary subroutines to initialize the problem space and 
scattering object(s) and for the incident waveform, far zone 
transformation, field update subroutines, outer radiation 
boundary conditions and field sampling. 

The main routine begins with the include statement and then 
appropriate data files are opened, and subroutines ZERO, BUILD 
and SETUP are called to initialize variables and/or arrays, build 
the object (s) and initialize the incident waveform and 
miscellaneous parameters, respectively. Subroutine SETFZ is 
called to intialize parameters for the near to far zone 
transformation if far zone fields are desired. 

The main loop is entered next, where all of the primary 
field computations and data saving takes place. During each time 
step cycle, the EXSFLD (TE) , EYSFLD (TE) , and EZSFLD (TM) 
subroutines are called to update the x, y, and z components of 
the scattered electric field. These scattered field equations 
are based upon the development given in [3]. RADEXY, RADEYX (TE) 
and RADEZX, RADEZY (TM) outer radiation boundary conditions are 
called next to absorb any outgoing scattered fields. Time is 
then advanced 1/2 time step according to the Yee algorithm and 
then the HXSFLD (TM) , HYSFLD (TM) , AND HZSFLD (TE) subroutines 
are called to update the x, y, and z components of scattered 
magnetic field. Time is then advanced another 1/2 step and then 
either near zone fields are sampled and written to disk in 
DATSAV , and/or the near zone to far zone vector potentials are 
updated in SAVFZ. The parameter NZFZ (described later) in the 
common file defines the type of output fields desired. 

After execution of all time steps in the main field update 
loop, subroutine FAROUT is called if far zone fields are desired 
to compute the far zone fields and write them to disk. At this 
point, the execution is complete. 


5 


SUBROUTINE SETFZ (TE, TM) 

This subroutine initializes the necessary parameters 
required for far zone field computations. The codes as furnished 
compute backscatter far zone field and can compute bistatic far 
zone fields for one scattering angle (i.e. one <p angle). Refer 
to reference [4] for a complete description of the two 
dimensional near to far zone transformation. Other versions of 
this subroutine provide for multiple bistatic angles. 

SUBROUTINE SAVFZ (TE, TM) 

This subroutine updates the near zone to far zone vector 
potentials . 

SUBROUTINE FAROUT (TE, TM) 

This subroutine changes the near zone to far zone vector 
potentials to far zone electric field 0 (TM) and <p (TE) 
components and writes them to disk. 

SUBROUTINE BUILD (TE, TM) * 

This subroutine "builds" the scattering object (s) by 
initializing the IDONE (TE) , IDTWO (TE) , and IDTHRE (TM) arrays. 
The IDONE- IDTHRE arrays are for specifying perfectly conducting 
and lossy dielectric materials. Refer to Figure 1 for a diagram 
of the basic two dimensional Yee cell for the TE and TM case. 

For example, setting an element of the IDONE array at some I,J 
location is actually locating dielectric material at a cell edge 
whose center location is I+0.5,J. Thus, materials with diagonal 
permittivity tensors can be modeled. The default material type 
for all ID??? arrays is 0, or free space. By initializing these 
arrays to values other than 0, the user is defining an object by 
determining what material types are present at each spatial 
location. Other material types available for IDONE-IDTHRE are 1 
for perfectly conducting objects and 2-9 for lossy non-magnetic 
dielectrics. It is assumed throughout the code that all 
dielectric materials are non-magnetic (i.e. the materials have a 
permeability of n 0 ) . This subroutine also has a section that 
checks the ID??? arrays to determine if legal material types have 
been defined throughout the problem space. The actual material 
parameters (e and a) are defined in subroutine SETUP. The 
default geometry is a 0.25 m radius, perfectly conducting 
cylinder . 

The user must be careful that his/her object created in the 
BUILD subroutine is properly formed. 

When it is important to place the object in the center of 
the problem space, NX and NY should be odd. This is due to the 



field locations in the Yee cell and also the placement of the E 
field absorbing boundary condition surfaces. 
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If the object being modeled has curved surfaces, edges, etc. 
that are at an angle to one or more of the coordinate axes, then 
that shape must be approximately modeled by lines and faces in a 
"stair-stepped” (or stair-cased) fashion. This stair-cased 
approximation introduces errors into computations at higher 
freguencies. Intuitively, the error becomes smaller as more 
cells are used to stair-case a particular object. 

SUBROUTINE DPLATE (TE) 

This subroutine builds squares of dielectric material by 
defining two each of IDONE and IDTWO components corresponding to 
one spatial square of dielectric material. It can also be used 
to define thin (i.e. up to one cell thick) dielectric or 
perfectly conducting wires. Refer to comments within DPLATE for 
a description of the arguments and usage of the subroutine. 

SUBROUTINE SETUP (TE, TM) * 

This subroutine initializes many of the constants required 
for incident field definition, field update equations, outer 
radiation boundary conditions and material parameters. The 
material parameters e and a are defined for each material type 
using the material arrays EPS and SIGMA respectively. The array 
EPS is used for the total permittivity and SIGMA is used for the 
electric conductivity. These arrays are initialized in SETUP to 
free space material parameters for all material types and then 
the user is required to modify these arrays for his/her 
scattering materials. Thus, for the lossy dielectric material 
type 2, the user must define EPS (2) and SIGMA (2 ) . The remainder 
of the subroutine computes constants used in field update 
equations and boundary conditions and writes the diagnostics 
file. 

SUBROUTINE EXSFLD (TE) 

This subroutine updates all x components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDONE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE EYSFLD (TE) 

This subroutine updates all y components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTWO array are used to determine the type of material present 
and the corresponding update equation to be used. 
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SUBROUTINE EZSFLD (TM) 

This subroutine updates all z components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTHRE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADEYX, RADEXY (TE) and RADEZX, RADEZY (TM) 

These subroutines apply the outer radiation boundary 
conditions to the scattered electric field on the outer 
boundaries of the problem space. 

SUBROUTINE HXSFLD (TM) 

This subroutine updates all x components of scattered 
magnetic field at each time step. The standard non-magnet ic 
update equation is used. 

SUBROUTINE HYSFLD (TM) 

This subroutine updates all y components of scattered 
magnetic field at each time step. The standard non-magnet ic 
update equation is used. 

SUBROUTINE HZSFLD (TE) 

This subroutine updates all z components of scattered 
magnetic field at each time step. The standard non-magnetic 
update equation is used. 

SUBROUTINE DATSAV (TE, TM) * 

This subroutine samples near zone scattered field quantities 
and saves them to disk. This subroutine is where the quantities 
to be sampled and their spatial locations are to be specified and 
is only called if near zone fields only are desired or if both 
near and far zone fields are desired. Total field quantities can 
also be sampled. See comments within the subroutine for 
specifying sampled scattered and/or total field quantities. When 
sampling magnetic fields, remember the £t/2 time difference 
between E and H when writing the fields to disk. Sections of 
code within this subroutine determine if the sampled quantities 
and the spatial locations have been properly defined. 

FUNCTIONS EXI, EYI (TE) and EZI (TM) 

These functions are called to compute the x, y and z 
components of incident electric field. The functional form of 
the incident field is contained in a separate function SOURCE. 


8 


FUNCTION SOURCE (TE, TM) * 

This function contains the functional form of the incident 
field. The code as furnished uses the Gaussian form of the 
incident field. An incident smooth cosine pulse is also 
available by uncommenting the required lines and commenting out 
the Gaussian pulse. Thus, this function need only be modified if 
the user changes the incident pulse from Gaussian to smooth 
cosine. A slight improvement in computing speed and 
vectorization may be achieved by moving this function inside each 
of the incident field functions EXI, EYI and so on. 

FUNCTIONS DEXI, DEYI (TE) and DEZI (TM) 

These functions are called to compute the x, y and z 
components of the time derivative of incident electric field. 

The functional form of the incident field is contained in a 
separate function DSRCE. 

FUNCTION DSRCE (TE, TM) * 

This function contains the functional form of the time 
derivative of the incident field. The code as furnished uses the 
time derivative of the Gaussian form of the incident field. A 
smooth cosine pulse time derivative is also available by 
uncommenting the required lines and commenting out the Gaussian 
pulse. Thus, the function need only be modified if the user 
changes from the Gaussian to smooth cosine pulse. Again, a 
slight improvement in computing speed and vectorization may be 
achieved by moving this function inside each of the time 
derivative incident field functions DEXI, DEYI and so on. 

SUBROUTINE ZERO (TE, TM) 

This subroutine initializes various arrays and variables to 

zero. 

V. INCLUDE FILE DESCRIPTION (TEACOM.FOR, TMACOM . FOR ) * 

The include files, TEACOM.FOR, TMACOM. FOR, contain all of 
the arrays and variables that are shared among the different 
subroutines. These files will require the most modifications 
when defining scattering problems. A description of the 
parameters that are normally modified follows. 

The parameters NX and NY specify the size of the problem 
space in cells in the x and y directions respectively. For 
problems where it is crucial to center the object within the 
problem space, then NX and NY should be odd. The parameter NTEST 
defines the number of near zone quantities to be sampled and NZFZ 
defines the field output format. Set NZFZ=0 for near zone fields 
only, NZFZ=1 for far zone fields only and NZFZ=2 for both near 
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and far zone fields. Parameter NSTOP defines the maximum number 
of time steps. DELX and DELY (in meters) define the cell size in 
the x and y directions respectively. The 0 incidence angle (in 
degrees) is defined by PHINC and the polarization is defined by 
the code that is being used. Parameters AMP and BETA define the 
maximum amplitude and the e temporal width of the incident 
pulse respectively. BETA automatically adjusts when the cell 
size is changed and normally should not be changed by the user. 
The far zone scattering angle is defined by PHIFZ. The codes as 
furnished perform backscatter computations, but this parameter 
could be modified for a bistatic computation. 

VI. SCATTERING WIDTH COMPUTATIONS 

Companion codes, SWTEA . FOR , SWTMA . FOR , have been included to 
compute scattering width versus frequency. Each code uses the 
file name of the FDTD far zone output data (FZOUTTE.DAT, 
FZOUTTM.DAT) and writes data files of far zone electric field 
versus time (FZTE.DAT, FZTM.DAT) and scattering width versus 
frequency (SWTE.DAT, SWTM.DAT). The scattering width 
computations are performed up to the 10 cell/A 0 frequency limit. 
Refer to comments within these codes for further details. 

VII. RESULTS 

As previously mentioned, the codes as furnished model an 
infinite, 0.25 m radius, perfectly conducting cylinder and 
compute backscatter far zone scattered field at an angles of 
0=180 degrees. 

Figures 2-4 show the TM far zone electric field versus time 
and the scattering width magnitude and phase for the 0.25 m 
radius perfectly conducting cylinder. Keep in mind that the fsr 
zone time domain electric field shown here is not the actual time 
domain scattered field. The actual far zone time domain 
scattered field can be obtained by an FFT of the FDTD time domain 
results, then multiplying the FFT by the appropriate frequency 
domain factor (described in [4]) and performing an inverse FFT. 

Figures 5-7 show the TE far zone electric field versus time 
and the scattering width magnitude and phase for the 0.25 m 
radius perfectly conducting cylinder. 

VIII. SAMPLE PROBLEM SETUP 

The codes as furnished model an infinite, 0.25 m radius, 
perfectly conducting cylinder and compute backscatter far zone 
scattered field at an angle of 0=180 degrees. The corresponding 
output data files are also provided, along with codes to compute 
scattering width using these data files. In order to change the 
code to a new problem, many different parameters need to be 
modified. A sample problem setup will now be discussed. 



10 


Suppose that the problem to be studied is scattering width 
backscatter versus frequency from a 0.30 m radius dielectric 
cylinder with a dielectric constant of 4e 0 using a 0-polarized 
field. The backscatter angle is 0=60.0 degrees and the frequency 
range is up to 3 Ghz . 

The incident field is 0-polarized which indicates the TM 
code must be used. Since the frequency range is up to 3 Ghz, the 
cell size must be chosen appropriately to resolve the field IN 
ANY MATERIAL at the highest frequency of interest. A general 
rule is that the cell size should be 1/10 of the wavelength at 
the highest frequency of interest. For difficult geometries, 

1/20 of a wavelength may be necessary. The free space wavelength 
at 3 GHz is A 0 =10 cm and the wavelength in the dielectric coating 
at 3 GHz is 5 cm. The cell size is chosen as 1 cm, which 
provides a resolution of 5 cells/A in the dielectric coating and 
10 cells/ A 0 in free space. Numerical studies have shown that 

choosing the cell size < 1/4 of the shortest wavelength in any 
material is the practical lower limit. Thus the cell size of 1 
cm is barely adequate. The cell size in the x and y directions 
is set in the common file through variables DELX and DELY. Next 
the problem space size must be large enough to accomodate the 
scattering object, plus at least a five cell boundary (10 cells 
is more appropriate) on every side of the object to allow for the 
far zone field integration surface. The default problem space 
size of 201 by 201 is adequate and provides a 75 cell border 
around the cylinder. As an initial estimate, allow 2048 time 
steps so that energy trapped within the dielectric layer will 
radiate. Thus parameter NSTOP is changed to 2048. If all 
transients have not been dissipated after 2048 time steps, then 
NSTOP will have to be increased. Truncating the time record 
before all transients have dissipated will corrupt frequency 
domain results. Parameter NZFZ must be equal to 1 since we are 
interested in far zone fields only. To build the object, simply 
change the RADIUS variable in the BUILD subroutine to 0.30. In 
the common file, the incidence angle PHINC has to be changed to 
60.0 respectively, and the cell sizes (DELX and DELY) are set to 
0.01. Since dielectric material 2 is being used for the 
dielectric coating, the constitutive parameters EPS (2) and 
SIGMA (2) are set to 4e 0 and 0.0 respectively, in subroutine 
SETUP. This completes the code modifications for the sample 
problem. 

IX. NEW PROBLEM CHECKLIST 

This checklist provides a quick reference to determine if 
all parameters have been defined properly for a given scattering 
problem. A reminder when defining quantities within the code: 
use MKS units and specify all angles in degrees. 
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TEACOM . FOR , TMACOM . FOR : 

1) Is the problem space sized correctly? (NX, NY) 

2) For near zone fields, is the number of sample points correct? 
(NTEST) 

3) Is parameter NZFZ defined correctly for desired field 
outputs? 

4) Is the number of time steps correct? (NSTOP) 

5) Are the cell dimensions (DELX, DELY) defined correctly? 

6) Is the incidence angle (PHINC) defined correctly? 

7) For other than backscatter far zone field computations, is 
the scattering angle set correctly? (PHIFZ) 

SUBROUTINE BUILD: 

1) Is the object completely and correctly specified? 

SUBROUTINE SETUP: 

1) Are the constitutive parameters for each material specified 
correctly? (EPS and SIGMA) 

FUNCTIONS SOURCE and DSRCE: 

1) If the Gaussian pulse is not desired, is it commented out and 
the smooth cosine pulse uncommented? 

SUBROUTINE DATSAV : 

1) For near zone fields, are the sampled field types and spatial 
locations correct for each sampling point? (NTYPE, IOBS, JOBS) 
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XI . FIGURE TITLES 

Fig. 1 Standard two dimensional Yee cell showing placement of 
electric and magnetic fields for the TE and TM case. 

Fig. 2 Far zone scattered field versus time for 0.25 m radius 
perfectly conducting cylinder using TM polarization. 

Fig. 3 Scattering width magnitude versus frequency for 0.25 m 
radius perfectly conducting cylinder using TM 
polarization. 

Fig. 4 Scattering width phase versus frequency for 0.25 m 
radius perfectly conducting cylinder using TM 
polarization. 

Fig. 5 Far zone scattered field versus time for 0.25 m radius 
perfectly conducting cylinder using TE polarization. 

Fig. 6 Scattering width magnitude versus frequency for 0.25 m 
radius perfectly conducting cylinder using TE 
polarization . 

Fig. 7 Scattering width phase versus frequency for 0.25 m 
radius perfectly conducting cylinder using TE 
polarization. 
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